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ABSTRACT 

In this paper we describe our convective hydrocodes for radial stellar pulsation. We adopt the 
KuhfuB (1986) model of convection, reformulated for the use in stellar pulsation hydrocodes. Phys- 
ical as well as numerical assumptions of the code are described in detail. Described tests show, that 
our models are numerically robust and reproduce basic observational constraints. 

We discuss the effects of different treatment of some quantities in other pulsation hydrocodes. 
Our most important finding concerns the treatment of the turbulent source function in convectively 
stable regions. In our code we allow for negative values of source function in convectively stable 
zones, which reflects negative buoyancy. However, some authors restrict the source term to non- 
negative values. We show that this assumption leads to very high turbulent energies in convectively 
stable regions. The effect looks like overshooting, but it is not, because turbulence is generated 
by pulsations. Also, turbulent elements do not carry kinetic nor thermal energy, into convectively 
stable layers. The range of this artificial overshooting (as we shall call it) is as large as 6 local 
pressure scale heights, leading to unphysical internal damping through the eddy-viscous forces, in 
deep, convectively stable parts of the star. 
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1. Introduction 

Pulsation hydrocodes play a key role in understanding the variability of clas- 
sical pulsators: 5 Cephei and RR Lyrae stars. They enable us to model the light 
and radial velocity curves and to study the modal selection in these stars. The first 
hydrocodes were purely radiative, as it was believed, that convection should not 
alter the pulsation properties of the stars, specially close to the blue edge. Indeed, 
radiative hydrocodes were very successful in reproducing many of the observed fea- 
tures (see e.g., Buchler 1998 and references therein). However several problems re- 
mained (Buchler 1998, Kovacs & Kanbur 1998), with modeling of the double-mode 
pulsations being the most severe among them. Success of convective hydrocodes in 
solving this longstanding problem (KoUath et al. 1998, Feuchtinger 1998) focused 
attention on modeling classical pulsators with convective hydrocodes. 
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Turbulent convection is an important physical process acting in many types of 
stars. It is essentially three dimensional process transporting energy through many 
length-scales: from macroscopic eddy cells to microscopic molecular scales, were 
energy is dissipated. Stellar convection acts in hard-turbulence regimes (extremely 
high Rayleigh numbers (> 10^^) and extremely small Prandtl numbers (~ 10~^, 
Gehmeyr & Winkler 1992a) which furthermore compUcates its modehng. How- 
ever, many essential features of convection may be described with simple one- 
dimensional models. Of these, the famous mixing-length theory (MLT) of Bohm- 
Vitense (1958) is most commonly used and underlies many other more complicated 
models. MLT however, is not suitable for stellar pulsation problems, since it is a 
local and time-independent theory. In pulsating variable stars, large-scale motions 
of gas interact with smaller-scale turbulent motions and time-dependent models are 
necessary to describe the coupling between them. Such one-dimensional models of 
turbulent convection were developed by many authors. Appropriate equations may 
be obtained through the Reynolds averaging technique (see e.g., Stanisic 1985). All 
quantities are decomposed into mean and fluctuating parts, for example for velocity 
and temperature we have U = U + U' , T = T + T' , respectively. Hydrodynamic 
equations decouple into mean and fluctuating (turbulent) equations, which are cou- 
pled by second order correlations, like U'U' or U'T' , that need to be modeled in 
order to close the system. Such procedure introduces several dimensionless, order 
of unity, free parameters, usually denoted by a-s. In one-equation models, turbu- 
lent equations are reduced to one equation for turbulent energy, = U'^/2. An 
excellent review showing in detail the derivations and approximations made in dif- 
ferent models, is given by Baker (1987). We will focus on two models that are 
commonly adopted in some recent linear as well as non-linear calculations, namely 
the Stellingwerf model (1982) and the models based on the work of KuhfuB (1986). 

Stellingwerf (1982) truncated the set of three turbulent equations derived by 
Castor (1968, unpublished) to a one-equation model for the turbulent energy, ap- 
plying MLT motivated closure relations. To model the small-scale turbulent dis- 
sipation, Stellingwerf introduced the eddy-viscous pressure term, in an ad hoc 
way. KuhfuB (1986) model is self-consistent, with all necessary modeling based 
on physical arguments. Eddy- viscous terms result from first-order modehng of the 
Reynolds tensor. All the correlations are modeled in a consistent way, using diffu- 
sion approximation. This leads to fully differentiable formulation, contrary to the 
Stellingwerf model, in which effects of buoyant deceleration of the turbulent eddies 
must be neglected in convectively stable regions. This leads to extreme overshoot- 
ing in Stellingwerf model as was shown by Gehmeyr & Winkler (1992a,b). These 
authors performed a detailed comparison of both models. They favour the KuhfuB 
model pointing other shortcomings of the Stellingwerf treatment. Some drawbacks 
of the original Stellingwerf theory are overcome by Bono & Stelhngwerf (1992, 
1994), who propose a better treatment of convectively stable regions (cf. Section 5). 
In our hydrocodes we adopt the KuhfuB model as it is self-consistent and based on 
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firm physical grounds. 

In this paper we present detailed description of our convective hydrocodes. In 
the first part of the paper (Sections 2 and 3) we discuss the turbulent convection 
model we adopt in our code: physical formulation in Section 2 and numerical im- 
plementation in Section 3 and Appendices. In the second part of the paper (Sec- 
tions 4 and 5) we present results of some tests of our hydrocode (Section 4), and we 
discuss the effects of different formulations used in other hydrocodes (Section 5). 
Conclusions are collected in Section 6. 



2. Physical description of the model 

Our convective hydrocodes implement the time-dependent turbulent convection 
model proposed, in its original form, by KuhfuB (1986). This model was reformu- 
lated by Wuchterl & Feuchtinger (1998) and Feuchtinger (1999) for the use in 
stellar pulsation hydrocode (Vienna code in the following). Also the hydrocode of 
Kollath et al. (2002) (Florida-Budapest code in the following) uses the model based 
on KuhfuB derivation (see however Section 5). We use a very similar formulation 
as in the above mentioned codes, however, we use some different assumptions and 
somewhat different parametrisation. Therefore, for clarity and completeness, we 
reproduce below all the equations and quantities we use in our codes. 

Momentum and energy equations are given by 

dU 1 3 GMr 

dE dV _ \ d[R\F,. + f::)] 

dt ^ dt~ p R^dR ' ^ ^ 

det^ pdV _ 1 djR^Ft 
'dt^ '~dt~~p R^dR 
Sum of the last two equations form the total energy equation 



+ Pt— = -- +£g + C. (3) 



In the above equations, U stands for the fluid velocity, which is time derivative of 
radius, R 

Mr is mass enclosed in radius /?, y is specific volume, which is inverse of density, 
p . E and P are pressure and energy of the gas including radiation, while and 
Pt are turbulent energy and turbulent pressure. Following Wuchterl & Feuchtinger 
(1998), we denote viscous energy and momentum transfer rates by Eq and Uq (note 
however, that Wuchterl & Feuchtinger 1998 use Uq = pUq and Eq = pEq). Fr, 
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Fc and are radiative, convective and turbulent flux, respectively. The C term 
describes the coupling between gas energy and turbulent energy equations. Below, 
we give detailed description of all model quantities. 

Turbulent pressure, , corresponds to trace part of the Reynolds tensor, while 
trace-free part leads to turbulent viscosity terms, Uq and Eq . These quantities are 
given by 



1 d 
4 1 fdU 



4nX 



_ 471 ax 

dM ' 



where 



and 



X = 



16 



d{U/R) 



is kinetic turbulent viscosity. A 

the pressure scale height. 

The couphng term, C, is given by 



= oHp is mixing-length, and Hp 



(6) 

(7) 

(8) 

(9) 
(10) 

-dR/d\nP is 



C = S-D-Dr 



(11) 



The D term describes the dissipation of turbulent kinetic energy into thermal en- 
ergy (turbulent cascade), while Dr describes the radiative cooUng of the convective 
elements (see Wuchterl & Feuchtinger 1998). Both these terms always damp the 
turbulent motions. The source function, S, is responsible for generation of turbu- 
lent energy through buoyant forces. It may drive as well as damp the convective 
motions. D, Dr and 5 are given by 



3/2 



(12) 



Dr=—^ —e„ 



TPQ 
cpHp' 



n. 



(13) 



(14) 



where a is Stefan-Boltzmann constant, cp is specific heat at constant pressure, 
Q = idV /dT)p is thermal expansion coefficient and K is opacity coefficient, n is 
correlation between entropy and velocity fluctuations, which is given by 



n = OJOLsel^^cpY. 



(15) 
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In the above expression, Y is superadiabatic gradient (dimensionless entrophy gra- 
dient), given by 

cp oR 

In convectively unstable layers, 7 > 0, part of the total energy is transported by 
convective flux, Fc 

F, = —pTU = aacpTcpey^. (17) 
a. 

The turbulent kinetic energy flux. Ft , is approximated by 

F, = -a,pAey'^. (18) 

In the adopted model it is the sole cause of the overshooting of turbulent eddies into 
convectively stable layers. 

As was shown by Wuchterl & Feuchtinger (1998), diffusion approximation 
used in modeling the IT correlation may be violated in some regions of the star. 
To fix this problem, these authors introduced the concept of convective enthalpy 
flux limiter, which we keep in our code as an option. In this case, expression (15) 
should be replaced by the following formula 



„ flw 1/2^ r f3T 

n = ^--e/F, sj-^-aa.cpY 



(19) 



where w = E + P /p is specific enthalpy, and F^ is the flux limiter function. 

For the radiative transfer one may use time-dependent, detailed treatment (Vi- 
enna code), or simpler models, like diffusion approximation (Florida-Budapest 
code). The first approach is certainly much more accurate and leads to better phys- 
ical description of the model's structure, specially in the outer parts. However, in 
classical pulsator's models, time dependent treatment gives essentially the same 
results as simple diffusion approximation (Kovacs & Kanbur 1998, Feuchtinger, 
Buchler & KoUath 2000). Thus, we adopt the diffusion approximation in our hy- 
drocode, which has the advantage of very low numerical costs. In this approxima- 
tion radiative flux, Fr, is given by 

F. = -^4./?^l^. (20) 

Radiation pressure, Pr, and radiation energy, Er, are included in P and E to- 
gether with gas contribution: P = Pg + Pr, E = Eg+Er. We have Pr = aT^/S and 
Er = aT'^/p, where a is radiation constant. Pressure, energy as well as other ther- 
modynamic quantities, are calculated as a function of T and V from the equation- 
of-state (EOS). We use either simple analytical EOS (Stelhngwerf 1982), or de- 
tailed EOS tables published by the OPAL group (Rogers et al. 1996). For the opac- 
ity coefficient we use the Rosseland mean. By default we use OPAL opacity tables 
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(Iglesias & Rogers 1996), generated for the solar mixture of Grevesse & Noels 

(1993) . At low temperatures these are supplemented by Alexander & Fergusson 

(1994) opacities (see Pamyatnykh 1999). 

Above formulae contain eight order of unity parameters. These are: a, a^, 
ttp , ttj , , ttc , ttd and Yr • Theory provides no guidance for their values, however, 
some standard values are in use. Parameters and (Xc were introduced by Yecko, 
KoUath & Buchler (1998), and are not present in the original KuhfuB (1986) deriva- 
tion, were = 2/3 and ttc = . Neglecting the radiative losses, local static so- 
lution of equations (l)-(3) may be reduced to MLT solution if ttc = = 1 /2-^2/3 
and a^; = 8/3-^2/3 (KuhfuB 1987). With radiative losses included in the model, 
exact MLT solution cannot be reproduced. We follow Wuchterl & Feuchtinger 
(1998) who opt for = 2a/3. We will refer to the above quoted values of a^, 
, ttc , ttd and Yr as standard values. We stress however, that they are not based 
on any firm physical considerations and therefore, should be treated as reference 
values only. In general, values of the a -parameters should be chosen to satisfy 
the observational constraints. However, results published up to date indicate, that 
no unique set of convective model parameters reproduces models of both RR Lyrae 
and Cepheids in stellar systems of different metallicities. In Table 1, we summarize 
the a parameters present in our model, and we give their relation to parametrisation 
used in Florida-Budapest code. 



Table 1 

Parameters of the discussed convection model. In ttie ttiird column we give a standard values, as 
described in Section 2. Fourth column gives the relation between our parametrization and 
parametrization used in Florida-Budapest code (barred alphas, KoUath et al. 2002). Only 
relation is not obvious, which results from complicated parametrisation of source function in 
Florida-Budapest code. No relation is given for the radiative losses, as this effect is treated 
differently in both codes (see Section 5.1) . For eddy viscosity treatment see Section 5.2. 



quantity a std. value relation to Florida- 

Budapest code 



mixing-length 


a 






a = 


= a 


eddy viscosity 












turbulent pressure 


ap 


2/3 


a,, 


= dp 


tiu-bulent source 




l/2i 


/273 


as 


= ds^Ud 


turbulent dissipation 




8/3i 






= 0.d 


convective flux 


a.c 


1/2, 




ac 




turbulent flux 


Of 






a. 


= cx, 


radiative losses 


Yr 


2y/3 
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3. Numerical representation of the model 

In pulsation hydrocodes, hydrodynamic equations are discretized on a mesli, 
wliicli is eitlier Lagrangean (fixed mass zones) or adaptive. More sopliisticated, 
adaptive grid is very useful in resolving the narrow features present in classical pul- 
sator's models such as shocks or hydrogen partial ionization regions (PIR). These 
are not very well resolved by Lagrangean mesh. However, the use of adaptive mesh 
in purely radiative models just smoothed the light curves, not changing their over- 
all shape (Buchler 1998). In convective models, hydrogen PIR is widened by the 
convective motions and therefore, is numerically less troublesome. Light curves 
are smooth already with the Lagrangean mesh. Finally, Feuchtinger et al. (2000) 
compared the results obtained with Florida-Budapest code (Lagrangean version) 
with the results obtained with Vienna code (adaptive mesh). Linear results as well 
as light/velocity curves agree very well. Therefore, instead of increasing numeri- 
cal costs by additional equation for adaptive mesh, we decided to keep the simple 
Lagrangean mesh. 

Our codes are based on radiative hydrocodes originally written by Stelling- 
werf (1975) with some later modifications (Kovacs & Buchler 1988). Inclusion 
of turbulent convection model however, requires significant changes in numerical 
methods that are used to solve the hydrodynamic equations. Below we describe the 
numerical schemes we use in our codes: in static model builder (Section 3.1), in 
the Unear code (Section 3.2) and in the nonlinear direct-time integration hydrocode 
(Section 3.3). By default, all these codes use the same analytical EOS (Stelhngwerf 
1982) and opacity procedures, the latter adopted from Warsaw-New Jersey stellar 
evolution code (Pamyatnykh 1999). 

3.1. Construction of the static model 

In classical pulsators, inner parts of the star do not participate in the oscilla- 
tions. We model the outer parts of the star, so-called envelope, only. We neglect 
rotation and magnetic fields. The model is specified by its mass, M, luminosity, L, 
effective temperature, T^q, and chemical composition, X and Z. We are not bound 
by evolutionary tracks in choosing these parameters. We also need to specify the 
a-parameters entering the convection model we use. 

The model is divided into mass zones. All quantities are defined either at 
the zones (thermodynamic quantities, T , V , P, £) or in between, at so-called 
interfaces (R, U , Y , fluxes; see Appendix A). Static model is constructed in two 
steps. In the first step, we construct an initial model without turbulent pressure and 
turbulent flux (a^ = = 0), with turbulent energy and coupling term defined at the 
interfaces. This initial model is constructed by integration of the static equilibrium 
equations from the surface inward. The final model, with turbulent pressure and 
turbulent flux included, and with turbulent energy and couphng term defined at the 
zones, is constructed through the multivariate Newton-Rhapson iterations. Below 
we briefly describe these two steps of model construction. 
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In the first step, we integrate tiie static, time-independent form of equations 
(l)-(4) {U = 0, d/dt = 0) from the surface inward. Initial conditions at the surface 
are: outer pressure, P = 0, radius of the outermost interface, Rn = {L/A'kgT^^YI'^ 
and the temperature of the outermost zone, = fT^^^. By default we use / = 
1/2 resulting from the Eddington approximation and keep f = '}>IA resulting from 
the exact solution of the gray atmosphere as an option. We decided to use the 
same mesh structure as was used in radiative codes (see e.g.,Kovacs & Buchler 
1988). Several outer zones, A'/i, have equal masses, DMN , down to the anchor 
zone, located at the hydrogen PIR. Temperature of this zone, T^, is fixed, Ta = 
1 lOOOK (or Ta = 15000 K, see Section 4). In the inner part of the envelope, below 
the anchor, masses of the remaining N — Na zones increase geometrically inward, 
with the common ratio h. Temperature of the innermost zone is Tin- Such mesh 
structure is obtained through adjustment of DMN and h during integrations. First, 
we integrate from the surface down to the anchor zone and adjust DMN to obtain 
desired temperature in the anchor zone. Then, we integrate from the anchor zone 
down to the inner boundary, h is adjusted to obtain the desired inner temperature. 

Tin. 

In the second step, results of integration serve as an initial guess for the multi- 
variate Newton-Rhapson iteration of the final model (Appendix A, eqs. (36)-(38)). 
Turbulent flux and turbulent pressure are turned on (if desired in the computed 
model) and turbulent energy and coupling term are redefined at the zones. We treat 
the iterations as converged, if relative corrections to temperatures, radii, and tur- 
bulent energies in all zones/interfaces are lower than 10"^*^ . To preserve the mesh 
structure, iterations are repeated several times with DMN and h being successively 
changed. DMN is adjusted to match the desired temperature in the anchor zone. 
Then, to assure the smooth transition from the upper part (zones of equal mass) to 
the lower part of the envelope (zone mass increasing geometrically inward), zone 
mass ratio h must be adjusted. It is done in such a way, that the mass of the en- 
velope below the anchor is not changed during the whole iterative procedure. As 
a result, the inner temperature of the envelope. Tin , changes, but only by several 
Kelvins. 

During the iterations we constrain the outer temperature to T^ = fL/ (47ia/?^) 
and allow for free adjustment of the outer radius, . This is fully compatible with 
the outer boundary condition for luminosity in the nonlinear code {cf. Section 3.3). 
Outer radius increases if turbulent pressure is turned on in the computed model. 

For some models with turbulent pressure included or with strong turbulent flux 
(ttf > 0.1 ) we encounter convergence difficulties. These are overcome by gradually 
increasing the and/or to the desired value during the iterations. 

In Appendix A, we give the numerical representation of equations solved through 
the Newton-Rhapson iterations and representation of all quantities that enter the 
static model computation. We also give representation for turbulent viscosity terms, 
that are not present in the static model, but enter the Unear stabiUty analysis and are 
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present in the nonlinear calculations. 
3.2. Linear stability analysis 

The original Stellingwerf (1975) codes implement the Castor (1971) numerical 
method in the linear analysis. As implementation of convection model into this 
scheme is very cumbersome, we drop this scheme and solve the full eigenvalue 
problem. 

We consider equations (l)-(3) and (5), however with energy equation (2) rewrit- 
ten in the following form 

dt V \dV)T) dt dM ' 

where cy is specific heat at constant volume. The static, not perturbed model is 
constructed by just described model buildefl We linearize this system, treating R , 
U , T and et as basic variables. We assume ~ exp(a?) dependence of the perturbed 
quantities, where a is the complex eigenvalue. Resulting eigenvalue problem is 
solved using canned eigenvalue solver as suggested by Glasner & Buchler (1993). 
We define the linear growth rates of the modes as the fractional growth of the ki- 
netic energy per pulsation period: r\ = 47i9?(a)/(0, where CO = 3(a) is the mode 
frequency. 

Additional equation for turbulent energy generates a new branch of eigen- 
modes. These are extremely damped with typical growth rates r\ < —I. Therefore, 
these modes are not expected to cause any troubles in nonlinear calculations (cf. 
Yecko et al. 1998). 

In addition to the frequencies and growth rates of the modes, linear code calcu- 
lates the radius, temperature, luminosity and turbulent energy eigenvectors. These 
are simply returned by the eigenvalue solver. Then, pressure work integrals may be 
simply calculated. For any pressure term, we have (e.g.. Castor 1971) 

rMR 

Wp{Mr) = -Ti 3{(5P)*(8V)}JM, (22) 

where integration is extended over the mass of the envelope. Also eddy viscosity 
contributes to the work integral. In Appendix C we derive following formulae for 
eddy-viscous work 

(■Mr ( 

Wev{Mr)=ti 3U5X) 

'However, we filter out the lowest turbulent energies as they generate numerical havoc, manifest- 
ing e.g., by erratic oscillations of the linear work integral in the inner envelope. We set e, = eg if 
e, <eQ, with eQ= i erg/g. We checked that results are independent of cq value at least in the range 
eo £ (10^^, 10^) erg/g. With ey below 10^^ erg/g numerical havoc appears. On the other side, setting 
the cutoff above 10^ erg/g we interfere with significant turbulent energies. 



6V 

/?3 



3^5^ 



dM. 



(23) 
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5X is easily calculated from already known eigenvectors (see eq. (9)). It is conve- 
nient to normalize the work integrals by the kinetic energy of the mode, Ek 



and such normalized work integrals are presented in the Figures of this paper. 
3.3. Nonlinear code 

In the nonlinear hydrocode, full set of nonlinear equations is integrated for- 
ward in time. The numerical scheme is very similar to the original radiative version 
(StelUngwerf 1975). Equations (1) and (3)-(5) are written in finite difference form, 

1/2 

using 05 = as the basic variable (Appendix B, eqs. (58)-(61)). In order to pre- 
serve the total energy of the envelope, time averages that appear in these equations, 
must be written in a careful way. We base our scheme (described in Appendix B), 
on the scheme of Fraley (1968) developed for radiative hydrocodes. 

Multivariate Newton-Rhapson iterations are used to solve the full nonlinear 
system of difference equations. Temperatures, 7] , radii, , and turbulent energies, 
et^i, at time t provide the initial guess, and values at time t + DT are iterated. 
We use constant time step, corresponding to roughly 600 steps per pulsation cycle. 
We treat the iterations as converged, if relative corrections to temperatures, bTi/Tf, 
and dU5i/U5i are smaller than 10^^ in all zones. Typically 3 to 6 iterations (up to 
^ 30 during the contraction phase) are necessary. For some models convergence 
difficulties are encountered. If convergence is not achieved in 60 iterations, the 
current iterations (not all computations) are restarted with halved time step. 

Difference equations are supplemented by the boundary conditions. At the 
inner boundary we have a rigid core of constant luminosity. Turbulent energy 
is set equal at the irmer boundary as well as in the outermost zone, et^ = 0. 
External pressure is set equal zero, and the outgoing luminosity, L^/, is given by 



From numerical point of view, turbulent energy equation requires special atten- 
tion. All its components depend on et in some power. Therefore, = is always a 
solution. Once et equals in some zone, it will stay equal 0, even if convective in- 
stability arises (see e.^., Gehmeyr & Winkler 1992a). Different numerical schemes 
were developed to overcome this problem (Yecko et al. 1998, Gehmeyr & Winkler 
1992b). We use the scheme very similar to that used by Yecko et al. (1998). Specif- 
ically, we add additional non-zero term in the turbulent energy equation, by slight 
modification of turbulent dissipation term, D (eq.(12)) 




(24) 



LN = 4nRlf-'aT^. 



3/2 



3/2 



D = ad 



A 



(25) 



where eo is a small, constant turbulent energy, which we set to eo = 10'*erg/g. In 
Section 5.3 we describe in detail how this correction works. 
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Inevitable component of radiative hydrocodes is artificial viscosity. It is nec- 
essary to handle shocks developing in the models. It acts as additional pressure, 
spreading the shock through several mass zones. As an unwanted by-product, ar- 
tificial viscosity limits the pulsation amplitude, being the main dissipative factor. 
Amplitudes of the radiative models depend on the parameters of the artificial vis- 
cosity, that should be adjusted to match the observational constraints. In convective 
hydrocodes, eddy viscosity provides physical source of dissipation. In principle, 
artificial viscosity is not necessary, however we keep it in our code. We use mod- 
ified Neumann-Richtmyer artificial viscosity (Stellingwerf 1975). It is described 
by two parameters, Cq which characterizes the strength of additional pressure, and 
by cut-off parameter, ttcut • Additional pressure turns on, only if relative speed of 
the consecutive zones exceed the local sound speed by fraction given by parameter 
ttcut ■ In our models we use Cq = 4 and high cut-off parameter, ttcut = 0.1. With 
such high cut-off, artificial viscosity does not turn on at all in most of the models. It 
always plays a subdued role, and is never present in the final hmit cycle pulsations. 

NonUnear hydrocode is supplemented with several data processing tools. For 
the converged limit cycle pulsations the nonlinear work integrals are calculated (see 
Appendix C). Also, bolomettic light and velocity curves are computed. During the 
pulsation cycle, photosphere sweeps through several Lagrangean zones. Therefore, 
photospheric values are extracted by interpolating to the exact black-body condi- 
tion, L = 4nR^oT'^ . Colour light curves are obtained through applying bolometric 
correction at each pulsation phase. We compute bolometric correction using static 
atmosphere models of Kurucz (2005). 

We stress that the nonhnear hydrocode is fully compatible with the linear one. 
Exactly the same Lagrangean mesh is used in both codes, as well as EOS and 
opacity procedures and numerical representation of all quantities. It is extremely 
important if one interprets nonhnear results in terms of linear ones. 

4. Tests of the code 

Some P Cephei models, computed with presented hydrocodes, were already 
published (Smolec & Moskalik 2007). In this Section we present hmited sample 
of other test calculations we have done. We stress that these are only test calcula- 
tions, not intended to model real stars. The goal is to show how our codes work, 
and to show that resulting models are reliable and numerically robust. We focus 
our attention on fundamental mode classical Cepheids. We consider models with 
two sets of convective parameters, A and B, given in Table 2. For a^, ttc and 
we use standard values (cf. Section 2). Set A represents the simplest convective 
model without turbulent pressure and turbulent flux, while in set B these effects are 
turned on. In both sets radiative losses are neglected as well as the turbulent flux 
hmiter. Again we stress, that we did not adjust convective parameters to match the 
observational constraints. For all the models discussed in this paper, we use Galac- 
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tic chemical composition, X = 0.7 and Z = 0.02, and following mass-luminosity 
relation: log(L/L0) = 3.56 log(M/MQ) + 0.7328 (Szabo etal. 2007). 

Table2 

Two sets of convective parameters considered in this work, a, , ttc , and are given in the 

units of standard values. 



Set 


a 




Us 




aa 






Yr 


A 


1.5 


0.20 


1.0 


1.0 


1.0 


0.0 


0.0 


0.0 


B 


1.5 


0.25 


1.0 


1.0 


1.0 


1.0 


0.01 


0.0 



Model envelopes are divided into 150 mass shells extending down to 2.5 • 10^ K. 
40 exterior zones have equal mass down to the anchor zone. Mass of the inte- 
rior zones increase geometrically inward. For models of set A we set the anchor 
temperature to llOOOK. Such anchor choice is not good for models of set B (and 
generally for models with turbulent flux) as the growth rates are not smooth along 
the sequence of models (e.g., for models of constant mass and differing tempera- 
ture). This is clearly numerical effect, resulting from poor resolution, as growth 
rates for models calculated with denser mesh (300 mass shells) exhibit smooth be- 
haviour. Smoothness may also be obtained by setting the anchor temperature to 
higher value. With 15000K growth rates are smooth. Therefore, we use this value 
of anchor temperature for models of set B. 

Static models constructed by model builder, are subject to linear stability analy- 
sis. This allows to determine the pulsation instability strips (IS) in the Hertzsprung- 
Russel diagram. For the discussed sets of convective parameters, IS are plotted in 
the left panel of Fig. 1 (set A) and in the right panel of Fig. 1 (set B). Thick solid 
lines define the fundamental mode IS, while thick dotted lines enclose first overtone 
IS. For set A, instability strips are shifted by ~200K toward higher temperatures, 
and are narrower than for set B. However, first overtone IS extends to slightly higher 
luminosities for set A. At low luminosities, the widths of the fundamental mode IS 
are roughly ~ 800K and ~ 950K for sets A and B, respectively. Fundamental mode 
IS widens toward higher luminosities for both sets of convective parameters. 

Structure of typical static model of set B is depicted in Fig. 2. This model has 
4.5M0 and lies inside the fundamental mode IS, 300K from its blue edge. Upper 
panel of Fig. 2 shows the run of superadiabatic gradient, 7 = V — Va, and adia- 
batic gradient, Va, versus the zone number. Arrows indicate the minima of Va 
connected with partial ionization regions (PIR). Hydrogen and both heUum PIR are 
clearly resolved. These regions give rise to convective instability, that is F > 0. 
In convectively unstable regions, part of the flux is carried by convection (solid 
Une in the lower panel of Fig. 2). In the discussed model turbulent flux (dotted 
Une in the lower panel of Fig. 2) diffuses the turbulent energies (dashed Une in the 
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lower panel of Fig. 2) also into convectively stable zones, where Y <0 (overshoot- 
ing). Displayed curves are smooth and all ionization features are properly resolved, 
which shows that our mesh is reasonable. 

Linear stability of the models may be studied through linear work integrals. For 
the discussed model, fundamental mode work integrals are plotted in Fig. 3. Total 
work is plotted with solid line, turbulent pressure work with dotted line and eddy- 
viscous work with dashed line. Upper panel of Fig. 3 shows the local work, that is 
work done in individual zones over one pulsation period, while lower panel shows 
the cumulative work integrals (expression (22)). Work integrals are normalized 
by kinetic energy of the mode. As a result, total cumulative work at the surface, 
is equal to the growth rate of the mode. Radiative damping acting in the interior 
of the model is overcome by the driving through the K -mechanism acting in the 
second helium and hydrogen-helium PIRs. Eddy viscosity has always stabilizing 
effect, while turbulent pressure work may contribute both to damping and driving. 
In the discussed model its overall effect is neutral (bottom panel of Fig. 3). 

Linear stability analysis provides information whether the model is stable or 
unstable against small perturbations. In case of instabiUty it tells nothing about 
the final pulsation state, its amplitude, light/velocity curves, or the final modal se- 
lection. This is the domain of nonlinear calculations. In these calculations static 
model is kicked with the scaled velocity eigenvector of the desired mode and time 
evolution of the model is followed. 

The consistency between linear and nonlinear calculations may be checked, 
by initializing nonlinear calculations in linear regime, that is with small surface 
velocity amplitude, say O.lkm/s. Then, the growth rate of the initialized mode may 
be calculated through the nonlinear work integrals (see Appendix C), as the average 
total envelope work, over several initial pulsation cycles. Such determination is not 
very accurate, as our initialization is never clean, that is in addition to the desired 
mode, also other modes, specially higher order, strongly damped overtones are 
present in the initial phase of integration. Nevertheless, such calculated growth 
rates (for both fundamental and first overtone modes) agree with the Mnear values, 
typically within ±3%, differences larger than ±5% appear only exceptionally. 

Full nonlinear calculations were performed for fundamental mode models, with 
convective parameters of set B. Integrations were carried over several hundred to 
few thousand pulsation cycles, till the Umit cycle pulsation was reached. In Fig. 4 
we show the nonlinear work integrals for the already discussed 4.5Mo model. 
These should be compared with linear work integrals displayed in Fig. 3. Artifi- 
cial viscosity does not contribute to the nonUnear work integrals at all. It is clearly 
visible in the lower panel of Fig. 4 that pulsation instability is saturated. Total cu- 
mulative work integral at the surface is equal to 0. In comparison to linear work 
integrals, sharp features are widened, being smeared by the motion of ionization 
fronts through the Lagrangean zones of the model. The strong damping by the 
eddy viscosity is clearly visible. In practice it means, that the parameter may 
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be used to control the limit cycle amplitude, similar to artificial viscosity parame- 
ters in radiative codes. Higher am , stronger the eddy- viscous damping and lower 
the amplitude. 

Fig. 5 shows the radial velocity (left panel) as well as bolometric light curves 
(right panel) calculated for sequence of fundamental mode models of set B. The 
Hertzsprung bump progression (e.g., Buchler, Moskalik & Kovacs 1990) connected 
with the 2:1 resonance between the fundamental mode and the damped second 
overtone, 2(0o = (O2 , is clearly visible. This resonance is crucial in shaping the 
light/velocity fundamental mode curves at periods around 10 days, were resonance 
center is located (Simon & Schmidt 1976, Kovacs & Buchler 1989). Considering 
velocity curves and starting from shorter periods (left panel of Fig. 5, bottom), the 
bump first appears at the descending branch of the velocity curve and then, as period 
of the fundamental mode grows (and P2/P0 ratio gets smaller), bump moves toward 
ascending branch and finally disappears. Our radial velocity curves are smooth, 
which is not the case for bolometric Ught curves during part of the expansion phase 
(Fig. 5, right). For massive models series of wiggles appears on the descending 
branch. These are due to poor resolution of our Lagrangean mesh. As envelope 
expands, very thin convection zone, sweeps through several mass shells, moving 
inward the model. 

Hertzsprung bump progression, as well as comparison of model curves with 
observed data, may be best studied through the Fourier decomposition parameters 
(Simon & Davies 1983). Here we focus on radial velocity curves. These are decom- 
posed into Fourier series, and amplitude ratios, Rki = A^/Ai , and Fourier phases, 
<t>)fei =^k — k^i > are calculated. In Fig. 6 we show the run of Ai , R21 and ^21 for 
computed radial velocity curves of set B, running 300K-600K to the red of the blue 
edge of the fundamental mode IS. Dots represent the observational data of Moska- 
Mk, Gorynya & Samus (2008, in preparation). Despite the fact that we haven't at- 
tempted to adjust the convective parameters to match the observational constraints, 
the overall agreement is quite good. Concerning Ai and /?2i , we note some prob- 
lems for shorter periods. Specifically, models do not reproduce the sharp increase 
of /?2i • Concerning ^21 , the model curves are shifted toward shorter periods in 
comparison to observations. This is coimected with the location of the 2cOo = CO2 
resonance center. The characteristic run of ^2\ is directly connected with this res- 
onance (Buchler, Moskalik & Kovacs 1990). Observationally the resonance center 
falls around 10 days. In numerical models, resonance location depends mostly on 
chosen mass-luminosity relation and on convective parameters. These were not 
adjusted to match the resonance center which explains the shift. In practice it is 
hard to infer the resonance location from observations only. Resonance location 
may be determined through the fit of the theoretical run of ^21 versus period to the 
observational data points (e.g., Kienzle et al. 1999). 
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5. Discussion of approximations and representations used in otiier 

hydrocodes 

In different hydrocodes, adopting essentially the same convection model, some 
processes or model quantities are modeled or treated in different ways. This con- 
cerns the modeling of radiative damping of the convective elements, modeling of 
the eddy-viscous terms and treatment of the turbulent source function and convec- 
tive flux in convectively stable regions. 

5.1. Modeling of radiative damping of the convective elements 

We describe the radiative damping of the convective elements through the ra- 
diative cooling term, as proposed by Wuchterl & Feuchtinger (1998) for the Vienna 
pulsation code (see Section 2). Different approach is adopted in Florida-Budapest 
code, where Peclet correction factor is used (Buchler & KoUath 2000). Peclet fac- 
tor multipUes the source term as well as convective flux, and accounts for decrease 
of convective efficiency in the limit of small Peclet numbers. As we have not im- 
plemented the Peclet correction in our hydrocode, we do not discuss the possible 
differences in computed models, caused by the two described treatments. 

5.2. Treatment of eddy viscosity 

In different pulsation hydrocodes, eddy-viscous terms are treated in different 
ways. We use the form derived by KuhfuB, resulting from first order modeling 
of the Reynolds tensor. The same form is used in the Vienna code (Wuchterl & 
Feuchtinger 1998). Many workers use eddy-viscous pressure, introduced in an 
ad-hoc way by Stellingwerf (1982). The form proposed by Kollath et al. (2002), 
slightly different from the original Stellingwerf (1982) form, is most commonly 
used (e.g., in Florida-Budapest code, Olivier & Wood 2005, Keller & Wood 2006). 
We will refer to these treatments of eddy viscosity as to KuhfuB eddy viscosity and 
Kollath eddy viscosity. If Kollath eddy viscosity is used, Eq and Uq terms are not 
present in equations (l)-(4), but additional pressure term of the form 

4 1/2 d /U\ 

should be placed next to turbulent pressure term, yielding following momentum 
equation 

dU 1 a , , GMr 

Simple algebra shows, that in the above equation the — term is dropped in 
comparison to momentum equation with KuhfuB eddy viscosity (eq.(l)). 

To check the effects of using eddy viscosity in the Kollath form, we imple- 
mented it in our hydrocode. Numerical representation of eddy-viscous pressure, 
Py , is given at the end of Appendix A. Numerical scheme is much simpler, than 
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in case of KuhfuB treatment, as Py is just additional pressure term, and the Fraley 
numerical scheme, the same as in purely radiative case, may be used (see Ap- 
pendix B). 

In Fig. 1 we compare instability strips for models with convective parameters of 
sets A and B (Table 2) computed with KuhfuB eddy viscosity (thick lines), and with 
KoUath eddy viscosity (thin lines). In the latter case instability strips are narrower. 
Comparing individual models we note higher growth rates for models computed 
with KuhfuB eddy viscosity. This effect may be compensated by lowering am in 
models computed with KoUath eddy viscosity (keeping other alphas fixed). 

There is no unique and representative way to compare the nonhnear models. We 
decided to compute the sequences of nonhnear models of set B, lying 300K from 
the blue edges of the instabihty strips, computed with KuhfuB and with KoUath 
eddy viscosities. As these linear edges do differ (Fig. 1) also effective temperatures 
of the computed nonlinear models of the same mass do differ. Fourier decomposi- 
tion parameters for the velocity curves are plotted in Fig. 7, sohd line for KuhfuB 
eddy viscosity, dotted hne for KoUath eddy viscosity. In agreement with hnear 
results, amplitudes are lower for models computed with Kollath eddy viscosity. 
Consequently, also /?2i is lower for these models. Concerning both treatments 
of eddy-viscosity give roughly the same results. Despite the lower amphtudes for 
models computed with Kollath eddy viscosity, quahtative run of Fourier parameters 
is the same in both treatments of eddy viscosity. 

5.3. Treatment of source function and convective flux in convectively stable 
regions 

As already mentioned by KuhfuB (1986), one of the shortcomings of the Stelling- 
werf (1982) theory is that the source function, 5, caimot damp the turbulent mo- 
tions, when a given layer becomes convectively stable during pulsation. This is 
due to the S ~ y/Y dependence, resulting from the chosen closing relation of 
the Stellingwerf's convection model (see Section 1). Similarly, for the convec- 
tive flux, Fc ~ \/y , in this theory. Such formulation leads to several problems, 
specially the range of overshooting is large and characteristic time scales for the 
growth and decay of turbulent energies cannot be defined in a reasonable fashion 
(Gehmeyr & Winkler 1992b). To overcome these problems. Bono & Stellingwerf 
(1992, 1994) modified the original Stellingwerf model and set 5 ~ sgn(y)-y/|y| 
and Fc ~ sgn(y)y^jF|. KuhfuB theory is void of these problems. It offers phys- 
ically well motivated, differentiable formulation, as we have 5 ~ 7 and ~ 7 
in this model. Thus, in convectively stable regions, both 5 and Fc have nega- 
tive values (see subsection 5.3.3 below). This is the formulation apphed in our 
code, and e.g., in Olivier & Wood (2005). However some of the workers using the 
KuhfuB model cut the source function, and in convectively stable regions, 7 < 0, 
set 5 = 0. The same restriction is apphed to convective flux. SymbohcaUy we 
write for such modified-KuhfuB model, 5 ~ 7+ and ~ 7+ . This is done e.g., \a 
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the Florida-Budapest hydrocode (KoUath et al. 2002). Buchler & KoUath (2000) 
show, that if a -parameters of the convective model are recaUbrated accordingly, 
StelUngwerf (1982) theory (S ~ \/y and Fc ~ \/Y) and modified-KuhfuB theory 
{S ~ y+ and Fc ~ 7+ ), as well as mixed formulation of Yecko et al. (1998) (5 ~ \/Y 
and Fc ^Y+) give very similar results. Below we argue that their conclusion re- 
sults from the fact, that in all these models source function is equal in convec- 
tively stable zones. If we do not limit the source function to non-negative values, 
both KuhfuB and S telling werf theories differ significantly, as studied analytically 
by Gehmeyr & Winkler (1992b). With our convective hydrocode modified to use 
non-negative source function we find further differences. Specially, we show that 
non-negative source function leads to significant turbulent energies in the inner, 
convectively stable zones. 

Structure of this Section is as follows. In subsection 5.3.1 we compare the 
results obtained with our code for two cases discussed above: 5 ~ 7 and ~ 7 
(our default formulation) and 5 ~ 7+ and Fc ~ 7+ (Florida-Budapest formulation). 
In subsection 5.3.2 we show, that in case of 5 ~ 7, treatment of convective flux 
plays only a minor role. In subsection 5.3.3 we summarize the physical arguments 
justifying treatment used in our code. 

5.3..1 5 ~ 7 and Fc ~ 7 versus 5 ~ 7+ and Fc ~ 7+ 

We focus our attention on the simplest models with convective parameters of set 
A. We discuss two convection models, namely our standard model with 5 ~ 7 and 
Fc ~ 7 , and modified one with 5 ~ 7+ and F^ ~ 7+ . We will refer to these con- 
vection models, as well as to pulsation models computed with them, as to NN and 
PP, respectively. It turns out that hnear results for both NN and PP models are al- 
most identical. Respective edges of the IS in the left panel of Fig. 1 would overlap, 
as they are shifted by less than 0.5K. Full nonlinear calculations were performed 
for sequence of models lying 300K to the red of the blue edge of the fundamental 
mode IS. We stress that respective models from NN and PP sequences have the 
same masses, luminosities and effective temperatures. Despite the fact that Unear 
results are almost identical, nonlinear results differ significantly. This is shown 
in Fig. 8, where we plot Fourier decomposition parameters, A\, R21 and (|)2i, for 
both NN (sohd Une) and PP (dotted line) models. The most striking difference 
concerns amphtudes, and consequently the /?2i ratio. For PP sequence amplitudes 
are significantly lower, specially at shorter periods. This leads to decrease of /?2i 
ratio, not observed for NN models. We trace these differences to very different be- 
haviour of the models in the deep, convectively stable (7 < 0) interior. To explain 
these differences we focus our attention on one model of 4.5M0 . Figs. 9 display 
the spatial profiles of turbulent energy during one pulsation cycle. In the upper 
panel (Figs. 9a,b) profiles of et for NN model are plotted, while in the lower panel 
(Figs. 9c,d) we plot the profiles for PP model. For each model two viewpoints are 
used to highUght the internal zones (panels a,c of Fig. 9) and external zones (panels 
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b,d of Fig. 9). Note different scales on vertical, logarithmic axes for NN and PP 
models. In both models internal zones are convectively stable. Y becomes positive 
in zones around 70, independently of the pulsation phase, and the model consid- 
ered (NN/PP). In Figs. 10 we show the corresponding nonlinear work integrals for 
NN model (upper panel of Fig. 10) and PP model (lower panel of Fig. 10). It is 
clear from Figs. 9 and 10 that in the PP model very high turbulent energies are 
present in convectively stable zones below zone 70 (zones 40-70), despite the fact 
that turbulent flux is not turned on. Such high turbulent energies cause significant 
eddy-viscous damping in the deep interior, as is well visible in the lower panel of 
Fig. 10. Consequently, amphtudes are lower for PP model. For NN model, tur- 
bulent energies are negUgible in the discussed internal zones (below zone 70), and 
hence the eddy-viscous damping is not present there (upper panel of Fig. 10). 

To explain the reasons for such differences we need to analyze the turbulent 
energy equation for NN and PP case in detail. We rewrite eq. (3) in the following 
form 



Turbulent flux and turbulent pressure terms are dropped (we consider set A of con- 
vective parameters), and we make use of eq. (25). The time derivative from the 
left-hand-side is to be balanced by the source term (5-term), turbulent dissipation 
term (D-term), correction term (eo-term), and eddy-viscous energy transfer term 
(£^-term). Eg and eg terms are always positive and thus, always drive the turbulent 
energies. D-term is always negative and thus, always contributes to the damping 
of turbulent energies. 5 -term may drive, as well as damp the turbulent energies in 
case of NN models, and always drive the turbulent energies in case of PP models. 

1/2 

It is also important to notice that 5 and Eq terms depend on et like ~ , while 

3 /2 

D-term depends on like ~ . We also note that eo-term has slightly differ- 

1/2 

ent form in Florida-Budapest code, where it is equal to (Xdeoe/ /A (Yecko et al. 
1998). However, the following discussion concerning PP models, does not depend 
on the exact form of eo -term, what we checked numerically (by using the eo -term 
in the form of Yecko et al. 1998) and what will become clear from the discussion 
below. 

We focus our attention on the bottom boundary of the convectively unstable 
region, where Y changes its sign (around zone 70), and on the internal, convectively 
stable regions. Below we show that the different treatment of 5-term in NN and PP 
convection models leads to the observed differences in computed pulsation models. 

In NN model, as Y becomes negative, also 5-term becomes negative and to- 
gether with D-term they damp the turbulent energies. As is visible in Fig. 9a, 
around zone 70, extremely rapid fall of turbulent energies, by roughly 25 orders of 
magnitude, happens. This effective damping is mainly due to the 5-term, and its 

1/2 

~ e/ dependence on et . D-term plays only a minor role in the described turbu- 




D— term 



eo^term 



(28) 
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lent energy fall. The driving effect of Eq -term is overcome by the damping effect 
of the 5 -term. Turbulent energies fall rapidly, and only the eo-term prevents them 
from falUng to zero. This is the only role of this term. We note in passing, that as 
5, D and Eq terms depend on et in some power, = would be a solution of 
equation (28). As et becomes in some zone, it will stay zero, even if convective 
instabiUty arises. Therefore, we need a non-zero, small term in turbulent energy 
equation, that would act as a seed for turbulent energies, as convective instability 
arises. In convectively stable regions, distribution of turbulent energies depends 
mainly on the balance between S, eo and Eq terms. D-term plays negligible role 

3 /2 

due to its ~ dependence on . As is well visible in Fig. 9a, with our choice of 
^0 (^0 = lO^erg/g), turbulent energies in the deep envelope (below zone 70), stabi- 
lize at a very small level, ~ 10^'^ — 10~^'*erg/g. This level depends on eo value, 
which should be chosen in order to assure very small, non-zero energies in con- 
vectively stable regions. This is the case with our choice. Distribution of turbulent 
energies at this low level is shaped by physics represented e.g., in the 5 -term, hence 
the small bump of around zone 40, caused by the iron opacity bump. Turbulent 
energies of order of 10^^^ — lO^'^^erg/g are negligible and physically not impor- 
tant. On the other side, they are sufficient to rebuilt the high turbulent energies if 
convective instability arises. This is clearly visible in Fig. 9b. As convective in- 
stability sweeps into external zones, turbulent energies are ignited from negligible 
level ( 10^'"erg/g) to full strength (lO'^erg/g). We stress again that the value of eo 
determines the very small level of turbulent energies in convectively stable zones 
of the model, but as such, it is not responsible for it. Without eo turbulent energies 
would fall to zero, which is, as mentioned, unwanted. In any case, eo -term does not 
fasten this fall of amplitudes as it is always a driving term. Therefore, a rapid fall 
of turbulent energies as Y becomes negative, is entirely caused by physical terms. 
Thus, eo-term plays only a numerical role, which we checked carefully by varying 
its value (in reasonable range). 

In the PP model situation is different. When Y becomes negative, source term 
is set equal to 0. The damping by the D-term reduces the turbulent energies. As is 
visible in Fig. 9c, turbulent energies around zone 70 fall roughly by three orders of 
magnitude from lO^^erg/g to roughly 10^'^erg/g. Then, the damping effect of the 
Z)-term is balanced by the driving effect of 5^ -term. Note that Z)-term depends 

3/2 1/2 

on et Uke ~ e/ , while Eg depends on et Uke ~ . As et falls, the damping 
strength of D-term falls more rapidly, than the driving strength of Eq -term. In 
the absence of damping 5 -term (which would depend on et in the same way as 
Eq) balance between D-term and fi^-term sets the turbulent energies on relatively 
high level, 10^ — lO^'^erg/g in the zones 40-70 (Fig. 9c). Then, a slow dechne 
of et below zone 40 is observed. This decline reflects the vanishing amplitude of 
the fundamental mode, as one moves inward the model. Consequently Eq-term 
decreases inward the model and vanish at the inner rigid boundary, where = eo 
(eq. (28)). eo-term plays negUgible role in the described balance and its value 
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(in reasonable ranges) has no effect on the described turbulent energy distribution. 
Even more. As we restart the nonUnear calculations for the full ampUtude model 
with eo -term dropped, we find no numerical problems, and no change in model 
behaviour is observed. In convectively stable zones, 40-70, turbulent energies are 
reduced by only 3 orders of magnitude compared to the values in the center of 
convective zone, and they are still high enough to produce significant eddy- viscous 
damping in the deep interior. This damping is clearly visible in the lower panel of 
Fig. 10, and is responsible for lower amplitudes of PP models visible in Fig. 8. 

It is clear from the presented discussion, that the different treatment of source 
term in convectively stable regions is responsible for the observed differences be- 
tween NN and PP models. The crucial point is, that in PP model 5 = in con- 
vectively stable zones. The form of source term in convectively unstable regions 
(5 ~ y+ vs. S '-^ is not very important as was checked by Buchler & KoUath 
(2000). Nonhnear work integrals presented by Buchler & Kollath (2000) display 
the same feature as visible in lower panel of Fig. 10 - strong eddy-viscous damping 
in more than 30, convectively stable model zones. If we allow for negative values of 
source function, strong eddy-viscous damping is not present in convectively stable 
zones. This result does not depend on the exact form of the source function, either. 
Both with 5 ~ 7 (our approach, upper panel of Fig. 10) and with 5 ~ sgn(7) 
(Bono & StelUngwerf 1994) eddy viscous damping is not present in the iimer, con- 
vectively stable parts of the model (see work integrals e.g., in Bono, Marconi & 
StelUngwerf 1999). 

As we described, high turbulent energies in the PP model, are partly caused 
by the turbulent energy driving through the Eg -term. This may be further con- 
firmed, by calculation of NN and PP models without eddy viscosity (a^ = 0). In 
the NN model situation in the convectively stable regions does not change. Tur- 
bulent energies are still extremely small. In the PP model, in convectively stable 
zones, 5 -term equals to 0, and et = eo is now a solution. With et = eo = lO^erg/g 
turbulent energies are 9 orders of magnitude smaller in comparison to turbulent 
energies in the center of convective zones (to be compared with 3 orders of mag- 
nitude reduction in case of / 0). Such energies are negligible. Eddy-viscous 
damping in the internal zones is not present. NN and PP models calculated with 
cCm = give qualitatively the same results. As eddy- viscous terms do not enter the 
static structure calculations, linear results are actually the same in both PP and NN 
models, as already mentioned. We have also checked that the presented discussion 
does not depend on the eddy viscosity form used in the model (see Section 5.2). If 
eddy-viscous pressure is used instead of Eg and Ug terms, quaUtatively the same 
turbulent energy profiles are observed in PP models. Turbulent energies are gener- 
ated through —Pv^ term which enters eq. (28) instead of ii^-term. 

In the above discussion we considered the models with convective parameters 
of set A, that is without turbulent pressure, and what is more important without 
turbulent flux, which is, in the adopted model, responsible for the overshooting. 
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Turbulent flux diffuses the turbulent energies beyond the convective instability re- 
gions, and hence significant turbulent energies may be present in convectively sta- 
ble zones. This is physical overshooting. In the just discussed models, turbulent 
flux was turned off. The NN models do not display significant overshooting, as 
turbulent energies are rapidly damped in convectively stable regions. Contrarily, 
PP models display high turbulent energies in convectively stable, internal zones, 
despite the fact that turbulent flux is turned off. This effect looks like overshoot- 
ing, but in fact it is different from overshooting and we will call this phenomenon 
artificial overshooting. Turbulent energies are generated at the cost of pulsations 
(through the fi^-term) and are not effectively damped, because of the neglect of 
buoyant forces. Also, turbulent elements do not carry kinetic energy (turbulent flux 
turned off) nor heat, as convective flux is equal to in convectively stable layers of 
PP model, by definition. Internal zones, with significant turbulent energies, lead- 
ing to the additional eddy-viscous damping, cover no less than 6-7 local pressure 
scale heights (zones 40-70 in the discussed models). Thus, the range of artificial 
overshooting is extremely large. 

Test calculations done with turbulent flux turned on show, that all the conclu- 
sions from the above discussion remain unchanged. This is supported by Figs. 1 1 
and 12. In Fig. 11 we show the Fourier decomposition parameters for models of 
set B. Solid fine is for NN models, while dotted fine for PP models. We observe 
the same differences as visible in Fig. 8 for models of set A. In Figs. 12 we plot 
the profiles of turbulent energy for 4.5M0 model of set B calculated with PP con- 
vection model. These are to be compared with Figs. 9c,d. Distribution of turbulent 
energies in the deep interior (Fig. 9c and 12a) is qualitatively the same for models 
with and without turbulent flux. It is still shaped by the balance of D and Eq -terms. 
The range of this artificial overshooting may be larger than in case without turbu- 
lent flux, but it cannot be smaller than just estimated value of ~(6-7) local pressure 
scale heights. Differences are visible in the external zones. Comparing Figs. 9d 
and 12b one sees smoother profiles of turbulent energies in model with turbulent 
flux, and higher turbulent energies in the outermost zones. These effects are caused 
by physical overshooting. 

There are also strictly numerical consequences of different treatments of the 
source function in NN and PP convective models. During the integration of PP 
model, we deal with few orders of magnitude in et , only (Fig. 9d). Convergence is 
relatively fast, and constant time-step may be used thorough the whole model inte- 
gration. Price to pay for the correct treatment of convectively stable regions in NN 
models, is greater numerical cost. NN models deal with turbulent energies spanning 
many orders of magnitude. If convective instability arises, turbulent energy must 
grow by several orders of magnitude, as is visible in Fig. 9b. This leads to conver- 
gence difficulties, and sometimes, the chosen time-step needs to be shortened, as 
discussed in Section 3.3. 
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5.3..2 Sr^Y and Fc^^Y versus 5 ~ F and F,. ~ 7+ 

From the discussion presented in the proceeding Section, it is clear that it is the dif- 
ferent treatment of the source function in convectively stable regions, that leads to 
the described differences in computed NN and PP models. Treatment of convective 
flux, that is Fc Y versus Fc F+ plays a minor role. This is easy to understand. 
PP and NN models differ in convectively stable regions. In PP model, convective 
flux is equal to zero in convectively stable zones by definition and in NN model, tur- 
bulent energies are very small there and thus, convective flux is neghgible anyway. 
This is fully supported by the model sequence calculated with 5 ~ 7 and Fc ~ }+ 
(model NP) - dashed line in Figs. 8 and 11. In Fig. 8 NP sequence is almost over- 
lapped with solid lines for NN sequence. Slightly higher difference between NN 
and NP models is visible in Fig. 11, where models with turbulent flux are displayed. 
Physical overshooting present in these models lead to higher turbulent energies in 
convectively stable zones, and hence effects of negative flux are stronger. However, 
the range of internal overshooting in NN and NP models is small, compared to PP 
models, and therefore, also the region with significant negative convective flux is 
very thin. We note that negative convective flux is always very small in our models, 
never exceeding 10 per cent of the total flux, usually being much smaller. 



5,3..3 Physical interpretation of negative 5 and 

It is clear from KuhfuB (1986) derivation, that there is no need to restrict the tur- 
bulent source function, as well as convective flux to non-negative values. Below 
we show, that source function is proportional to forces acting on turbulent eddy 
during its motion. As these forces act both in convectively stable, and convectively 
unstable zones, there is no justification to neglect them in the former case. 

Our reasoning follows the MLT considerations, concerning the acceleration of 
convective elements. For details of the derivations, we refer the reader to Cox & 
Giuli (1968, §14.3). Below we assume that the convective eddy moves adiabati- 
cally, which corresponds to the neglect of radiative losses. 

For the turbulent eddy, equation of motion may be written as follows 

r = -g--— = -g—, (29) 
p or p 

where Ap is excess density of the turbulent element. Thus, acceleration of the 
element, resulting from buoyant forces is equal to 

Ap 

a = -gy- (30) 

Excess density may be expressed through the excess temperature, AT , (Ap /p) = 
—TpQ{AT /T) which, after traveling the mixing length A, is 

Ar(A) = -A^^(V-Va). (31) 
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Using above relations and equation of hydrostatic equilibrium to eliminate g from 
eq. (30), we obtain a final equation for eddy acceleration 

This is a buoyant acceleration of turbulent eddy displaced up by distance A. It 
is positive in convectively unstable region (7 > 0) and negative in convectively 
stable region (7 < 0). Comparison with definition of source function (14), leads to 
following relation 

5 = a/J^a ~ Y. (33) 

Thus, the source function is proportional to the acceleration of convective eddies 
caused by buoyant forces. This acceleration does not vanish in the convectively 
stable regions (it becomes negative). Therefore, there is no reason to set the source 
function to zero in convectively stable regions of the stellar model. To the contrary, 
source function has to be negative, because buoyant forces are now slowing down 
motion of the convective eddies. Assuming 5 = in convectively stable regions 
amounts to neglecting buoyancy, which is physically incorrect. 

Overshooted elements carry kinetic energy (turbulent flux), as well as thermal 
energy (convective flux), which are dissipated in convectively stable zones. It is 
easy to show, that in convectively stable regions of the star, convective flux has to 
be negative, just as the source function. Indeed, when convective eddy overshoots 
down from the envelope convective zone, it becomes hotter than the surrounding 
medium (because Va > V). When eddy overshoots up, it becomes cooler than 
the surrounding medium. In both cases, convective flux is directed downwards 
{Fc < 0). 

For further support for negative source function and negative turbulent flux we 
refer the reader to KuhfuB (1986), Gehmeyr & Winkler (1992a,b), Canuto (1997). 

6. Conclusions 

In this paper we describe, our convective hydrocodes for radial stellar pulsation. 
Convection model we use, is based on the KuhfuiS (1986) model. In the first part 
of the paper we briefly describe the model and list all model equations and neces- 
sary quantities. Technical details, concerning numerical schemes and methods we 
use, are fully described in Section 3 and in Appendices. Many tests we have done, 
part of which are briefly described in Section 4, prove that our code works prop- 
erly. Computed models are numerically robust and reproduce well basic features 
observed in classical Cepheids, Uke Hertzsprung bump progression. 

There are several other hydrocodes, that adopt similar convection model as 
our code. However, in different hydrocodes different treatments of some quan- 
tities, such as the turbulent source function, S, or eddy- viscous terms are used. 
Consequences of these differences were not fully studied up to date. In Section 5 
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we compare some of these treatments. Our most important finding, concerns the 
treatment of source function in convectively stable zones. In our hydrocode, we 
allow for negative source function in convectively stable zones, which reflects neg- 
ative buoyancy, and is physically well motivated. However, in some other codes 
(e.g., Florida-Budapest code), source function is restricted to non-negative values. 
This corresponds to the neglect of buoyant forces in convectively stable layers of 
the model. Similarly, convective flux is also restricted to non-negative values. We 
find that such approach has several serious drawbacks, that we list below: 

(/) Due to neglect of negative buoyancy effects (by assuming non-negative 
source function), significant turbulent energies are present in convectively stable 
layers of the model. When negative source function is absent in convectively stable 
layers, balance between eddy-viscous driving and turbulent dissipation sets the tur- 
bulent energies at relatively high level, 10^ — 10^"erg/g. We call this phenomenon 
artificial overshooting, as turbulent energies are generated by pulsations. Also, tur- 
bulent eddies do not transport heat into convectively stable layers, as convective 
flux is equal to 0. 

(//) The range of this artificial overshooting is very large. Significant turbulent 
energies extend to more than 6 local pressure scale heights below the envelope 
convection zone. 

(Hi) Significant turbulent energies in the deep, convectively stable parts of the 
model, lead to strong eddy-viscous damping, and consequently to lower pulsation 
amplitudes of the models in comparison to models computed with negative source 
function. 

(iv) Physical overshooting, due to turbulent flux, plays a minor role, and de- 
scribed, high turbulent energies are present in convectively stable layers, regardless 
if turbulent flux is included in the model or not. 

In the next paper (Smolec & Moskalik 2008) we will show further conse- 
quences of assuming non-negative source function. We show crucial role of this 
assumption in double-mode Cepheid models computed with the Florida-Budapest 
code (KoUath et al. 2002). We will explain in detail the mechanism that leads to 
double-mode behaviour and we will show, that this mechanism does not work if 
source term is treated properly, that is, if we allow for negative values of the source 
function. 
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7. APPENDIX A: Numerical representation of the model: static case 

The model envelope is divided into mass zones, separated by the interfaces. 
In our notation both zones and interfaces are denoted by integers. Adjacent inter- 
faces of zone / have indices i — 1 (bottom interface) and i (upper interface), just as 
in the scheme below: 



j+1 

i+l 

/ M, DM2, U, R, Hp, Y, L^, L„ L„ 

i DM, T, V, Q, CP, K, P, P„ Eg, S, D, Dr 
i-\ 

Outermost zone, as well as outermost interface, have index A^. All quantities are 
defined either at the zones or at the interfaces, just as presented in the scheme 
above. For some quantities spatial averages, denoted by curly brackets, need to 
be calculated. Average of the zone quantity (e.g., temperature, T) is defined at 
the interface and average of the interface quantity {e.g., radius, /?) is defined at the 
zone, Uke in examples below 

{T]i = Q.5{Ti + Ti+,), 

{R}i = 0.5{Ri.,+Ri). 

We also need to calculate the spatial differences, which are numerical represen- 
tation of derivatives. We denote them by A. Spatial difference of zone quantity 
(e.g., pressure, P) is defined at the interface, and spatial difference of the interface 
quantity (e.g., luminosity, L) is defined at the zone, Uke in examples below: 

APi = Pi+i-Pi, (34) 



ALi = Li-Li_i. (35) 

Model builder solves the static version of equations (l)-(4). These are rewritten 
in a Lagrangean form. Instead of internal energy equation (2), we solve the total 
energy equation (4). In the static limit this equation reduces to total luminosity 
conservation condition, L = const . Complete set of equations is as follows 

.APi + APfi GM, 

= -AnRj '7,^ ^ , (36) 

' DM2i /?? ' 

* I 

= ^ + ^ + ^-\, (37) 

±^ 
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0=-5i^+C,. (38) 

Turbulent energy equation is defined at the zones, while luminosity conservation 
and momentum equations at the interfaces. Mass enclosed by radius Ri is denoted 
by Mi , mass of the zone i by DM,- and mass associated with the interface i by 
DM2,. We have 

DM2,- = 0.5 (DM,- + DM,-+i ) . (39) 

For given temperature, 7] , and volume, Vi , EOS procedure calculates other ther- 
modynamic quantities defined at the zones: pressure, P, , and energy, £, (both con- 
taining gas and radiation contribution), specific heat at constant pressure, cpj , and 
thermal expansion coefficient, (2,. Opacity procedure calculates the opacity, k, , 
also defined at the zone. Remaining zone quantities are turbulent pressure (eq. (6)) 

Pt,i = ^, (40) 
and the quantities entering the coupling term, C, (eqs. (12)-(14)) 

„ TiPiQij n \ 1/2 

3/2 

Dri = ^^et ,-. (43) 

Other quantities are defined at the interfaces. Pressure scale height is 

Numerical representation of the superadiabatic gradient (eq. (16)) is based on the 
formula given by Stellingwerf (1982) 

''. = ^#((^l(A/'.)-(>°S^«-log^,)Y (45) 



DM2,- {y},- VI CP. 
Numerical representation of the n term is 

n,- = aas{cp}iYi, (46) 

or if the flux Umiter is turned on 

ni = !FiFL ^ , (47) 
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C48, 



where: 



3 [ T 

Qi = aas{cp}iYi. (49) 

1/2 

Note, that in comparison to equations (15) or (19) we dropped the term , which 
now appears separately in the definitions of source term, Si , and in the convective 
luminosity definition below, where it is averaged at the interface. Our motivation 
was to assure that all quantities entering the coupling term in zone / , C; , depend on 
the turbulent energy in this zone only, . For the luminosities we have 

L,, = 4;r/?^^|^}n,W/^}, (50) 

C , N 3/2 _ 3/2 



A, = --aa,(4<)^//.,|^|^^:^g^, (51) 

^ 4a (47r/??)^7;.%/K,+i-7;.VK,- 
3 DM2; 1 _ iog(K,^fi/Ki) 

Averaging scheme in the expression for radiative luminosity comes from StelUng- 

werf (1975). 

Eddy- viscous terms are not present in the construction of static envelope, since 
they depend on U . However, they appear in the linear as well as nonlinear code. For 
completeness we present their numerical representation below. JJq term is defined 
at the interface, while Eq at the zone 



Eq i = 4nXi , (53) 



Ui/Ri-Ui-i[Ri-i 
DMi 

jnXi^ (54) 
where X, is zone quantity defined as follows 

Xi = ^;raa.i,.;f (55) 

As an option alternative to Uq and Eq terms, we implemented the eddy- viscous 
pressure, defined, Uke other pressure terms, at the zone 

16 1/2 1 , , 1 (Ui Ui-A 

P\i = - -rrT^^^me, '■ {R Hp}i — — — - - — . (56) 
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8. APPENDIX B: Numerical representation of tlie model: nonlinear scheme 

Let Z denote any physical quantity entering our model. In our notation Z^f^ 
stands for the value of Z, at some particular moment of time, denoted by upper 
index (n). After the time step DT , its value is Z^''^'\ Time difference will be 

denoted by capital D, for example Z){7, = uj"^^^ — J//"^ . We also need to calculate 
the average value of Z, during the time step, which we denote by (Z,) . Usually we 
set 

<Z,-)=^J"+V(1-^)ZW, (57) 

with ^ = 1 corresponding to fully implicit treatment and ^ = corresponding to 
fully explicit treatment. For some quantities, more complicated time averages are 
necessary, as required by energy conservation (see below). 

Finite difference form of equations (1), (4), (3) and (5) is then 

£^^,n{Ri}^M±P^ + GMjL)-(U„}^0, (58) 
DT \ DM2,- \Rj / 



DT 

D{Ei + Ui}) + {{Pi) + {Pt,i))DVi + ^A((L„-) + (L,,,-) + (L,,» -DT{E,,i) = 0, 

(59) 

DT 

Dm] + {Pt,i)DVi + ^ A(Lf,,-) - £)r(£^,,-) - DT{Q) = 0, (60) 

Rf'^'^=Rf'^+DT{Ui). (61) 

Numerical scheme used in nonlinear calculations, must preserve the total en- 
ergy during the integrations. Appropriate scheme for the radiative hydrocodes was 
proposed by Fraley (1968). We consider the total energy of the model envelope, 
equal to 



Etot = 



DMi{Ei + 05f ) + 0.5DM2iUf 



(62) 



and its change during one time-step, DEjqj /DT . Energy conservation requires 
DEtctt/DT = 0. Neglecting the eddy-viscous terms, Eq and Uq, convective equa- 
tions may be reduced to the form exactly corresponding to purely radiative case, 
with the following substitutions: E = E , for the total energy, P = P + Pj, 
for the total pressure and L = Lr + Lc + Lt, for the total luminosity. Thus, without 
eddy-viscous terms, following time-averagings are necessary in order to preserve 
the total energy (Fraley 1968) 

{Ui) = lut^'^+^-ut\ (63) 
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,2(n+l)^^(«+l)^(«)^^2(„)^^ (64) 

.^)^^n)^^- ^^^^ 
' i i 

The way we average pressures and luminosities has no effect on energy con- 
servation. In our code we adopt the values used in radiative codes, that is we set 
9 = 9, = 1/2 and Wr = Wc = Wt = 2/3 in the following 

{Pi) = ep^"+^^ + {i-e)p^"\ (66) 
(^M> = e.^!r^' + (i-e0^y, (67) 

<L,,-)=w4fi' + (l-w,)Ly, (68) 

{L,j) = w,L^:;'^ + {\-w,)Li% (69) 

<l,,)=w,l;;;+^)+(i-w.)l!J. (70) 

Taking eddy-viscous terms into account and adopting averaging scheme as de- 
scribed above, one finds, that the total energy changes during one time-step, by 



DEtot ^ 



DMi{E,,i) + 0.5DM2i{U,,i)2{Ui) 



(71) 



Thus, total energy is conserved (the above sum vanishes), if following definitions 
for (^Uq) and (S^) are adopted 

, , 4ji (Xi+i)-(Xi) 
fe)^4.(.,> <^')/W-^<^-)/<'"-) . (73) 

(x,) = e,x,'"+'i + (i-e.)x<">, (74) 

<R,) = P«!"+" + (l-P)«S">. (75) 

9„ and [3 are not restricted by energy conservation. We set 9[, = 1 and P = 1/2. 

Turbulent energy equation (60) is decoupled from the energy conservation anal- 
ysis. For the coupling term we write 

<C,)=^("+i) + (l-Y)cf"\ (76) 

and set Y = 1 • Values of 9„ and y were chosen experimentally to assure fast con- 
vergence of the nonUnear iterations. 



30 



A. A. 



9. APPENDIX C: Eddy-viscous work integrals 

In this Appendix we present the derivations of non-linear and linear eddy- 
viscous work integrals. We follow the derivation method presented in Unno et al. 
(1989). Following these authors we rewrite momentum equation (1) in the vector 
form. For eddy- viscous terms we use relations (7) and (8). Momentum and energy 
equations are 

^(^ + P„^ = _ivf+4«», (78) 
at dt p oR 

where ir is unit vector in the radial direction, (|) is gravitational potential, Ur is 
radial (only) component of the velocity vector U . Ptot = P + Pt F = Fr + 
Fc + Ft . We multiply the momentum equation by pf/ to obtain the equation of 
conservation of mechanical energy 

pM = _£iV.„-p£iV* + 4.p^«. (79, 
dt R am 

Adding this equation to energy equation multiplied by p, and using continuity 
equation {p{dV/dt) = VU) we get 

Integrating over the mass of the envelope we get 

where Etot is equal to the total energy of the envelope 

rM 



-tot 

'Mo 



pM 

/ {E + et + U^/2 + ^/2)dM. (82) 

JMn 



Note that the last integral vanish as at the inner boundary Ur = and at the outer 
boundary X = 0. The remaining integral is rewritten in the following form 

dE 1 r 

_J21 = _/ _^(p)dM- PtotUds, (83) 

dt JMo P Jr=Ro 

where Rq is radius of the star. The last surface integral vanish due to our outer 
boundary condition (Ptot = 0)- Using energy equation again, we obtain 



^£tot _ r'' 

dt JMo 



d{E + et) dV .^d{UR/R) 
+ "tot— -4JtA- 



dt '"'dt dR 



dM. (84) 
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Nonlinear work integral is obtained through integrating over the whole pulsation 
cycle, which yields 



W 



dt 







M 



dM 



Mo 



dV 

Ptot—-4nX 
dt 



d{UR/R) 
dR 



dt 



M 



dM 



Mo 



dV 

J tnt — En 



dt 



In the linear approximation we have 

rP rM 



W 



f dt dM 

Jo JMo 



dt 



dR 



(85) 



(86) 



Using relation, Ur = d{dR)/dt, and assuming 6z = 9t(5ze""') for the perturbed 
quantities, after laborious but straightforward algebra we arrive at 



pM rM 

W = -n 3 [(SPtot)* (8V)] dM + n 3 

JMn JMo 



{5xy 



/SV _3VdR 



dM. (87) 



The first integral correspond to ordinary pressure work integral, and the second 
term is eddy-viscous work integral, if KuhfuB form of eddy viscosity is used. In 
case of Kollath form, eddy viscosity has a form of ordinary pressure, and is simply 
included in Ptot • 
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log r.„ log 

Fig. 1 . Linear instability strips calculated for models with convectivc parameters of set A (left panel) 
and set B (right panel). Solid lines limit the fundamental mode instability strip, dotted lines enclose 
the first overtone instability strip. Thick lines refer to models calculated with our standard eddy 
viscosity form, while thin lines refer to models calculated with KoUath eddy viscosity form, discussed 
in Section 5.2. 
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Fig. 2. Static structure for 4.5M0 model of set B, lying 300K to the red of the fundamental mode 
blue edge. All quantities are plotted versus the zone number. Surface at right. In the upper panel 
we display the run of 7 = V — Vg and Vg . Labeled arrows mark the minima of V;, connected with 
partial ionization regions of indicated element. In the lower panel relative convective and turbulent 
fluxes are plotted, together with scaled turbulent energy. 
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Fig. 3. Linear work integrals versus the zone number, for model of Fig. 2. Local work in the upper 
panel and cumulative work in the lower. 
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Fig. 4. Nonlinear work integrals versus the zone number, for model of Fig. 2. Local work in the 
upper panel and cumulative work in the lower. 
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Fig. 5. Full amplitude radial velocity curves (left panel) and bolometric light curves (right panel) for 
models of set B, running parallel to the blue edge of the fundamental mode IS, 300K from it. Model 
masses are increasing by O.SM© , starting from 4.0M© at the bottom of the figures up to 9Mq at 
the top. Consecutive radial velocity curves are shifted by 25 km/s to allow comparison. Curves are 
labeled by Unear fundamental mode periods, Pq , and linear P2/P0 period ratios. 
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Fig. 6. Fourier decomposition parameters of radial velocity curves for models of set B, running par- 
allel to the blue edge of the fundamental mode IS, 300K, 400K, 500K and 600K from it. Amplitudes 
are scaled by constant projection factor equal to 1.4. Individual curves for sequence running 300K 
apart from blue edge are presented in Fig. 5. Dots represent observational data (Moskalik, Gorynya 
& Samus 2008, in preparation). 
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Fig. 7. Fourier decomposition parameters of radial velocity curves for models of set B, running 
parallel to the blue edge of the fundamental mode IS, 300K from it. Amplitudes are scaled by constant 
projection factor equal 1.4. Solid line for models with KuhfuB eddy viscosity, dotted line for models 
with Kollath eddy viscosity (see Section 5.2). 
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Fig. 8. Fourier decomposition parameters of radial velocity curves for models of set A, running 
parallel to the blue edge of the fundamental mode IS, 300K from it. Amplitudes are scaled by constant 
projection factor equal 1.4. Solid line for NN models, dotted line for PP models, and dashed line for 
NP models (see Section 5.3). 
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Fig. 9. Profiles of turbulent energy during one pulsation cycle for models discussed in Section 5.3. 
Panels a) and b) for NN model highlighting the internal and external parts of the model, respectively, 
while panels c) and d) present corresponding profiles for PP model. 
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Fig. 10. Nonlinear, local work integrals versus the zone number, for 4.5M0 model of set A, lying 
300K to the red of the fundamental mode blue edge. Upper panel for NN model, lower panel for PP 
model. 
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Fig. 11. Fourier decomposition parameters of radial velocity curves for models of set B, running 
parallel to the blue edge of the fundamental mode IS, 300K from it. Amplitudes are scaled by constant 
projection factor equal 1.4. Solid line for NN models, dotted line for PP models, and dashed Une for 
NP models (see Section 5.3). 




Fig. 12. Profiles of turbulent energy during one pulsation cycle for PP model with turbulent flux 
turned on, discussed in Section 5.3. 



